%% STA code for checking connections in paired recording, should work with whole cell or on cell spikes.
clear all
% dataFolder='10_03_17_01';


recNumb=[28];% 13 14 15 16]; %recording number.5
preCh=1; %input presynaptic channel.


if preCh == 0
    postCh=1;
else
    postCh=0;
end

% path1='C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data\2017\';
% path2=[ dataFolder(1:2) dataFolder(6:8) '\' dataFolder(1:8) '\' dataFolder '\'];
%path1='C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data\2018\08_18\08_24_18\';
 path1=['C:\Users\Tom�s\Desktop\Temp\'];
%path1=['C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data\2018\08_18\08_16_18\'];

postTrace=[];
preTrace=[];

for i=1:numel(recNumb)
    
    dataStruc(1)= IBWread([path1  'ch' num2str(preCh) '_' num2str(recNumb(i)) '.ibw']);
    dataStruc(2)= IBWread([path1 'ch' num2str(postCh) '_' num2str(recNumb(i)) '.ibw']);
    
    Dx=dataStruc(1).dx;  %take any IBW to find sampling rate.
    Fs=1/Dx;
    
    
    preTrace=[preTrace ; dataStruc(1).y]; %Horizontally concatenate data vectors
    postTrace=[postTrace; dataStruc(2).y]; %same as above
    
    
end


rawTraces(:,1)= preTrace;
rawTraces(:,2)=medfilt1(postTrace,20);



figure;plot(dataStruc(1).y(1:1e5))


prompt = 'set threshold for pk detection         ';
thresh = input(prompt)
close gcf


%% extract spike indices using peakfind
if thresh>0
[pks,locs]=findpeaks(rawTraces(:,1),'minPeakprominence',thresh,'minpeakDistance',500);
else
[pks,locs]=findpeaks(rawTraces(:,1)*-1,'minPeakprominence',-thresh,'minpeakDistance',500);
end
    
spikeIndex=locs(10:end-10);
window=50*1e-3/Dx;% first term is window in ms
for i=1:numel(spikeIndex)
    mat4STAPost(i,:)=rawTraces([(spikeIndex(i)-window):(spikeIndex(i)+window)],2);
    mat4STAPre(i,:)= rawTraces([(spikeIndex(i)-window):(spikeIndex(i)+window)],1);
end


%remove trials with more than 1 presynaptic spike (>10 Hz)
% 
% multipleSpikeTrials=0;
% for i=1:numel(spikeIndex)
%     pks2=[];
%     locs2=[];
%     [pks2 locs2]=findpeaks(mat4STAPre(i,:),'minPeakprominence',thresh,'minpeakdistance',200);
%     if numel(locs2)>1
%         mat4STAPre(i,:)=NaN;
%         mat4STAPost(i,:)=NaN;
%         multipleSpikeTrials=multipleSpikeTrials+1;
%     end
% end
% 
% multipleSpikeTrials


STATracePost=nanmean(mat4STAPost,1);
STATracePre=nanmean(mat4STAPre,1);
time=(1:numel(STATracePost))*Dx*1e3;
STATracePostMedFilt=medfilt1(STATracePost, 10);

fig2=figure(2);
fig2.Renderer='Painters';

subplot(2,1,1);

% plot(time(1e2:end),medfilt1(mat4STAPost([1:500],1e2:end),10),'color',[0 0 0 0.05], 'linewidth',1)
% hold on;
plot(time(1e2:end),STATracePost(1e2:end),'color',[ 130 0 159]./255,'linewidth',3)

ylabel('I (pA)');
set(gca,'fontsize',12)
set(gcf,'color','white');
set(gca,'box','off')
xlim([40 80])
xticklabels([])
title(recNumb)

subplot(2,1,2);

plot(time(1e2:end),medfilt1(mat4STAPre([1:500],1e2:end),10),'color',[0 0 0 0.01], 'linewidth',1)
hold on
plot(time(1e2:end),STATracePre(1e2:end),'color',[0 191 255]./255,'linewidth',3);

ylabel('Vm (mV)');
xlabel('Time (s)');
set(gca,'fontsize',12)
set(gcf,'color','white');
set(gca,'box','off')
xlim([40 80])

%% analysis of paired recordings 
%First analyze the size of the STA using 3000 spikes.
STABaseline=nanmean(STATracePost(600:950));
[STAPeak delay]=max(STATracePost(1000:1100));
delay=delay*Dx*1000;
STASize=STAPeak-STABaseline;
% 
% %Next analyze the frequency of IPSC and PC firing rate
% %use recording length and number of PC spikes to find mean firing rate
% 
meanPreFiringRate=numel(locs)/(numel(dataStruc(1).y)*Dx);
% 
% %Find IPSC frequency.
prom= 5
[pks2,locs2]=findpeaks(rawTraces(:,2),'minPeakprominence',prom,'minpeakDistance',50);

IPSCRate=numel(locs2)/(numel(dataStruc(1).y)*Dx);

%suptitle([path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Pair ' path2(end-2:end-1)]);

%% analyze success vs failure using the first derivative

succ=[]
for i=1:numel(spikeIndex)
    tempTrace=mat4STAPost(i,:);
    tempdiff=diff(medfilt1(tempTrace,50));
    if nanmean(tempdiff(1020:1050))>0.2
        succ(i)=1;
    else
        succ(i)=0;
    end
end


successTrials=mat4STAPost(logical(succ),:);
failureTrials=mat4STAPost(~(logical(succ)),:);

% determine amplitude of events 
for i=1:size(successTrials,1)
    successAmplitude(i)=max(successTrials(i,1050:1060))-mean(successTrials(i,990:1010));
end

for i=1:size(failureTrials,1)
    failureAmplitude(i)=max(failureTrials(i,1050:1060))-mean(failureTrials(i,990:1010));
end

% make histogram of success and failure amplitudes
figure
hold on
histogram(failureAmplitude,'BinLimits',[-10 50],'BinWidth',3,'faceColor',[0 0 0])
histogram(successAmplitude,'BinLimits',[-10 50],'BinWidth',3)
temp=gca
Ytick2=temp.YLim(2);
xlim([-10 50])
set(gca, 'YTick', [0:Ytick2/3:Ytick2 ]);
yticklabels({0,'','',Ytick2})
set(gca, 'XTick', [0 25 50])
set(gca, 'TickDir', 'out')
xlabel('amplitude (pA)')
ylabel('number of events')

mat4STAPostBS(i,:)=mat4STAPost(i,:)-prctile(mat4STAPost(i,:),10);

%run this section to check if classifier works well
% 
% for i=1:sum(succ)
%     figure
%      hold on
%     plot(nanmean(successTrials,1),'r','linewidth',2)
%     plot(successTrials(i,:),'k','linewidth',2)
%     vline(1000,'--r')
%     hline(mean(successTrials(i,990:1010)))
%     hline(max(successTrials(i,1050:1060)))
%     pause
%     close gcf
% end
% 
% for i=1:(numel(succ)-sum(succ))
%     figure
%     hold on
%     plot(nanmean(successTrials,1),'r','linewidth',2)
%     plot(failureTrials(i,:),'k','linewidth',2)
%     vline(1000,'--r')
%  hline(mean(failureTrials(i,990:1010)),'--b')
%     hline(max(failureTrials(i,1050:1060)),'--g')
%     pause
%     close gcf
% end





% for plotting success trial overlay
figure;plot(time,nanmean(successTrials,1))
figure;hold on
for i=1:sum(succ)
    plot(time,successTrials(i,:),'color',[0 0 0 0.1])
end
plot(time,nanmean(successTrials,1),'r','linewidth',2)
set(gcf,'color','white');
axis off
axis tight
xlim([40 70])
obj=scalebar('XLen',2,'XUnit','ms','YLen',20,'YUnit','pA','Position', [60 140])
title(['success trials = ' num2str(size(successTrials,1)) '/' num2str(numel(spikeIndex))]);
successRate=sum(succ)/numel(succ)
% 
figure;
hold on;
clear i
for i=1:size(failureTrials,1)
    plot(time,failureTrials(i,:),'color',[0 0 0 0.1])
end
    plot(time,nanmean(failureTrials,1),'r','linewidth',2)
set(gcf,'color','white');
axis off
axis tight
xlim([40 70])
obj=scalebar('XLen',2,'XUnit','ms','YLen',20,'YUnit','pA','Position', [60 140])




%% plot all postsynaptic traces for STA in semitransparent black
% % 
% figure;
% hold on;
% for i=1:numel(spikeIndex)
%     plot(time,mat4STAPost(i,:),'color', [0 0 0 0.05])
% end
% xlim([40 80])
% % 
% pkprom=5
% clear amplitude location
% for i=1:numel(spikeIndex)
%     clear pks locs 
%    [pks,locs]=findpeaks((mat4STAPost(i,:)-mean(mat4STAPost(i,:)))*-1,'minPeakprominence',pkprom,'minpeakDistance',20,'minpeakwidth',20);
%    if numel(pks)==0  % check if there are any IPSCS in trace
%      amplitude(i)=NaN; 
%      location(i)=NaN;
%    else if any(locs>1000 & locs<1100)
%     logicalIdx=locs>1000 & locs<1100; %select only pks that fall after the presynaptic spike and before the decay of the average ipsc.
%     amplitude(i)=max(pks(logicalIdx)); %if more than one peak exists, pick the largest one.
%     location(i)=max(locs(pks==amplitude(i))); %pick the corresponding location for idx.
%     
%        else 
%      amplitude(i)=NaN; 
%      location(i)=NaN;
%        end
%    
%    
%    
%   
%     end
%     
%    [amp(i) idx(i)]=max((mat4STAPost(i,(1000:1100))-mean(mat4STAPost(i,:)))*-1);
%      
% end

figure;
hold on
for i=1:numel(spikeIndex)
   
plot(time,(mat4STAPost(i,:)),'color',[0 0 0 0.1])
mat4STAPostBS(i,:)=mat4STAPost(i,:)-prctile(mat4STAPost(i,:),10);
% plot(mat4STAPost(i,:)-prctile(mat4STAPost(i,:),10),'color',[0 0 0 0.1])
end
plot(time,nanmean(mat4STAPost,1),'color',rgb('deepskyblue'),'linewidth',2)

% plot(nanmean(mat4STAPostBS,1),'color',rgb('deepskyblue'),'linewidth',2)


% figure;
% hold on;
% for i=1:numel(spikeIndex)
%     if ~isnan(amplitude(i))
%         plot (time,mat4STAPost(i,:)-mean(mat4STAPost(i,:)),'color',[0 0 0 0.1])
%     else
%        
%     end
% end

% % 
% figure;
% hold on;
% for i=1:numel(spikeIndex)
%     if isnan(amplitude(i))
%         plot (time,mat4STAPost(i,:)-mean(mat4STAPost(i,:)),'color',[0 0 0 0.1])
%     else
%        
%     end
% end


% set(gca,'fontsize',12)
% set(gcf,'color','white');
% set(gca,'box','off')
%   
    
% 
% figure; plot (time,mat4STAPost((~isnan(amplitude)),:)) %this works for some reason.
% 
% figure; hist((location(~isnan(amplitude))-1000)*1000*Dx,200);
% % 
% figure;
% shadedErrorBar((time),mat4STAPost((~isnan(amplitude)),:),{@mean,@std},'lineprops','-k','transparent

fclose all


%% New section to determine latency

% STApeakdelay
STAPeakDelay=(find(STATracePost==max(STATracePost(1000:1100)))-1000)*Dx*1e3
derivative=diff(STATracePost);
thresh=(5*std(derivative(900:1000)))+nanmean(derivative(900:1000));

%find index of derivative where it crosses the threshold after 1000 (1000
%is where the peak of the action potential happens).

test=derivative>thresh;
test2=test(1000:1200);
meanLatencySam = find(test2, 1, 'first')
meanLatencyMs=meanLatencySam*Dx*1e3;

meanLatencyMs


D = finddelay(STATracePre,STATracePost)
for i=1:sum(succ)
    delay(i)= finddelay(STATracePost,successTrials(i,:)); % Thgis makes no sense! it's just repeating the same thing over and over. GG
    
end

delayMs=delay.*Dx*1e3;
%this is ugly code but I need to get rid of those 0 value trials, so i'm
%going to just delete them.
delayMsNonZero=delayMs(delayMs>0);
figure;histogram(delayMs,'binwidth',0.1)

figure;hold on
for i=1:(numel(succ)-sum(succ))
%     amplitudeFail(i)=mean(failuremat4STA(i,1045:1050));
    plot(failuremat4STA(i,:),'color',[0 0 0 0.05])
end
%     
successMat4STA=mat4STAPostBS(logical(succ),:);  
failuremat4STA=mat4STAPostBS(~logical(succ),:);
figure;hold on
for i=1:sum(succ)
   plot(successMat4STA(i,:),'color',[0 0 0 0.1])
    amplitudeSucc(i)=mean(successMat4STA(i,1045:1055));
end



